home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YARON.C < prev    next >
C/C++ Source or Header  |  1988-11-28  |  2KB  |  90 lines
  1. /********************** YDblVdp for double vanderPol *************************/
  2. /********************* (C) 1986,7,8 by JAMES A. YORKE ************************/
  3.  
  4. #include "yinclud.h"
  5.  
  6. aron(k, Y)        /* 2 coupled vanderpol oscillators -- see also map 23 */
  7. double  k[],
  8.         Y[];
  9. {
  10.     double  bb,
  11.             tt,
  12.             t,
  13.             X,
  14.             W,
  15.             yy,
  16.             N,
  17.             exp(), expY, expXYW, expX, expWY, expN;
  18.  
  19.     t = Y[0];
  20.     X = Y[1];
  21.     W = Y[2];
  22.     yy = Y[3];
  23.     N = Y[4];
  24.     bb = C2;
  25.     tt = 1.+ beta * cos(C4 * t);
  26.     expY = -(rho - 1.) * C1 * tt * exp(yy);
  27.     expXYW = (C3 + C1) * tt * exp(X + yy - W);
  28.     expX = C1 * rho * bb * exp(-X);
  29.     expWY = (C5 + C1) * exp(W - yy);
  30.     expN = C1 * bb * exp(-N);
  31.     k[0] = 1.;
  32.     k[1] = expY + expX - C1;
  33.     k[2] = expXYW - (C3 + C1);
  34.     k[3] = expWY - (C5 + C1);
  35.     k[4] = expN - C1;
  36. #ifdef GARBAGE
  37.     if(num_lyap > 0) {
  38.         for(i = 0; i < num_lyap; i++) {
  39.             base = lyapzero + vec_dim * i;
  40.             k[base] = Y[base + 1];
  41.             k[base + 1] = () * Y[base]
  42.                 + C1 * (1.- X1 * X1) * Y[base + 1]
  43.                 + rho * Y[base + 2];
  44.             k[base + 2] = Y[base + 3];
  45.             k[base + 3] =
  46.                 rho * Y[base]
  47.                 + (-2 * C5 * X2 * Y2 - C6 - 3 * C7 * X2 * X2) * Y[base + 2]
  48.                 + C5 * (1 - X2 * X2) * Y[base + 3];
  49.         }
  50.     }
  51. #endif
  52. }
  53.  
  54. init_aron() {
  55.     extern int      aron();
  56.  
  57.     zeroth = 1;
  58.     vec_dim = 3;        /* the dimension of the Lyapunov vectors =
  59.                    phase space dim; the component y[4] is not
  60.                    being counted in this because it drifts
  61.                    slowly */
  62.     num_lyap = 0;        /* the number of Lyapunov numbers to be
  63.                    computed <= vec_dim */
  64.     lyapzero = 5;        /* y[lyapzero] is the zeroth coord of the
  65.                    zeroth lyapunov vector */
  66.     dim = lyapzero + num_lyap * vec_dim;
  67.  /* needed for rungekutta */
  68.     X_upper =.3;        /* x scale */
  69.     X_lower = -.3;
  70.     Y_upper = 1.5;        /* y scale */
  71.     Y_lower = -5.;
  72.     X_coord = 1;
  73.     Y_coord = 3;
  74.     C1 =.02;
  75.     C2 =.6667;
  76.     C3 = 35.84;
  77.     C4 = twopi;
  78.     C5 = 100.;
  79.     beta =.1;
  80.     rho = 18.;
  81.     steps_per_cycle = 200;
  82.     step = (twopi / C4) / steps_per_cycle;
  83.     DE = aron;        /*  DE is a pointer to a function   */
  84.     y[eqn0 + 1] = -0.08960;
  85.     y[eqn0 + 2] = -1.90847;
  86.     y[eqn0 + 3] = -1.90918;
  87.     y[eqn0 + 4] = -0.00662;
  88.     setequal(0, 1, eqn1);    /* reinitializes y[] */
  89. }
  90.